Biomolecular interactions determine how transcription factors recognize their DNA binding sites, how proteins interact with each other, and consequently how a biological system functions. Since many biological molecules bear considerable electric charge, electrostatic interactions are among the most important when studying biomolecular interactions. However, electrostatic interactions in biological systems are difficult to calculate accurately in practice. Aside from the significant charges carried by biomolecules such as DNA and proteins, the solvent itself, namely, water produces considerable electrostatic effects. Furthermore, hydrogen bonds, known to be involved in helix formation in both DNA and proteins, are essentially electrostatic in origin. Indeed, it seems that electrostatic effects often drive the physical-chemical processes in biological systems and, thereby, determine biological function. Therefore, any attempt to perform molecular dynamics (MD) simulations of biological systems will require an adequate description of these electrostatic forces. Previously, we have developed a rigorous method to calculate the crucial electrostatic forces in a biomolecular system. Exact analytical results have been obtained for systems with sufficient symmetry. For example we have completely analyzed a system consisting of a pair of objects with planar surfaces containing embedded charges, which is an idealized model of molecular contact. We have also developed methods to compute the electrostatic forces for a biomolecular system in which the atoms are represented by spheres. The method is rigorous in the context of the model and the accuracy can be tuned to any desired level. However, the underlying idea we developed to compute the electrostatic forces in a biomolecular system can, with some modification, be applied to a more realistic and flexible model in which a biomolecule is represented by an arbitrary surface which is decomposed into a set of small triangular patches. There is, of course, fundamental effect that require more than electrostatics to understand. For example, the short distance electronic interactions call for quantum mechanical understanding. Recently, we have developed a quantum mechanical approach that is suitable for many electron systems by finding the density functional needed for the density functional theory. To see how we can integrate our new approach with existing ones, we have also try to use currently available density functionals to perform calculations for peptide fragmentation probability. The experience seems encouraging and the computational results were written and published. For the past year, the main effort, however, was to implement the frameworks we developed in computer codes which, after optimization, will be used to demonstrate the utility of the method in a biomolecular system. The electrostatic code, which takes into account only classical effects, has finally converged into a workable prototype. Our next goal is to speed up the code to make it useful in molecular dynamics simulation by providing the correctly calculated molecular forces. Of course, to take into quantum mechanical effects, we have also been coding so as to use the newly devised universal density functional to compute the energy of the homogeneous electron gas, a long-standing computational problem using Feynman diagrams or similar techniques. The crux lies in the numerical integration of high dimensions, where the integrand may change sign frequently thus hinders the numerical accuracy required. Specifically, in terms of classical level code, we are currently testing the accuracy and adequacy of surface triagularization by comparing the exact analytical result with those numerically obtained via surface triangularization. To better explain the core idea of our surface charge method in a more approachable way, we have written a new article, recently published in American Journal of Physics, to provide intuitive understanding.